stanEx1
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}stan basic
normal distribution
ex1.stan
normal distribution
data{
int N;
vector[N] y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~normal(m,s);
}mdl=cmdstan_model('./ex1.stan')
y=rnorm(10,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.42 1.79 2.29 2.17 3.16 3.86par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y') 
data=list(N=length(y),y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -178.954
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 -7.14868 0.00640488 0.000400752 1 1 17
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.3 seconds.mle## variable estimate
## lp__ -7.15
## m 2.17
## s 1.24mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -8.05 -7.71 1.12 0.83 -10.35 -6.97 1.00 746 687
## m 2.17 2.18 0.52 0.47 1.30 2.99 1.00 647 609
## s 1.55 1.47 0.46 0.38 0.99 2.35 1.00 1154 1077mcmc$metadata()$stan_variables## [1] "lp__" "m" "s"mcmc$metadata()$model_params## [1] "lp__" "m" "s"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -8.05 -7.71 1.12 0.83 -10.35 -6.97 1.00 746 687
## m 2.17 2.18 0.52 0.47 1.30 2.99 1.00 647 609
## s 1.55 1.47 0.46 0.38 0.99 2.35 1.00 1154 1077
poisson distribution
y=rpois(20,3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 2.00 4.00 3.40 4.25 8.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=length(y),y=y)ex2-1.stan
in case fit poisson dist. to normal dist.
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> m;
real<lower=0> s;
}
model{
y~normal(m,s);
}
generated quantities{
vector[N] y1;
for(i in 1:N)
y1[i]=normal_rng(m,s);
}# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -29.1211
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -22.0597 0.00121336 0.000370207 1 1 10
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -22.06
## m 3.40
## s 1.83
## y1[1] 3.25
## y1[2] 0.83
## y1[3] 4.25
## y1[4] 5.06
## y1[5] 1.06
## y1[6] 3.24
## y1[7] 5.54
## y1[8] 1.40
## y1[9] 5.38
## y1[10] 6.61
## y1[11] 1.67
## y1[12] 4.02
## y1[13] 2.19
## y1[14] 7.61
## y1[15] 5.49
## y1[16] 4.65
## y1[17] 6.51
## y1[18] 4.31
## y1[19] 5.76
## y1[20] 1.75mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -21.28 -20.98 1.04 0.75 -23.47 -20.27 1.01 704 1032
## m 3.41 3.40 0.45 0.43 2.70 4.13 1.00 1628 1102
## s 2.02 1.98 0.37 0.35 1.52 2.68 1.00 796 1104
## y1[1] 3.29 3.29 2.16 2.05 -0.25 6.80 1.00 1889 1794
## y1[2] 3.32 3.32 2.14 2.11 -0.22 6.80 1.00 1879 1817
## y1[3] 3.41 3.37 2.10 2.00 0.03 6.89 1.00 1983 1931
## y1[4] 3.37 3.39 2.12 2.03 -0.18 6.93 1.00 2204 2012
## y1[5] 3.43 3.47 2.06 1.98 0.08 6.73 1.00 1699 1787
## y1[6] 3.48 3.53 2.10 1.99 -0.01 6.82 1.00 1951 1956
## y1[7] 3.36 3.42 2.10 1.97 -0.22 6.79 1.00 2156 1949
## y1[8] 3.37 3.35 2.12 2.00 -0.11 6.90 1.00 1963 1816
## y1[9] 3.35 3.37 2.06 1.96 -0.05 6.70 1.00 1978 1646
## y1[10] 3.41 3.34 2.13 2.01 0.06 6.85 1.00 1925 1975
## y1[11] 3.51 3.48 2.14 2.02 0.06 7.06 1.00 1798 1532
## y1[12] 3.38 3.37 2.11 2.06 -0.11 6.90 1.00 2274 2035
## y1[13] 3.44 3.47 2.11 2.03 -0.09 6.98 1.00 2130 1933
## y1[14] 3.38 3.37 2.07 2.07 -0.06 6.58 1.00 1823 1712
## y1[15] 3.47 3.48 2.15 2.09 -0.19 6.96 1.00 2006 1840
## y1[16] 3.42 3.46 2.12 1.99 -0.04 6.82 1.00 1887 1858
## y1[17] 3.42 3.40 2.10 1.95 -0.04 6.81 1.00 1896 1894
## y1[18] 3.41 3.44 2.12 2.04 -0.13 6.79 1.00 1867 1741
## y1[19] 3.44 3.39 2.05 2.01 0.21 6.85 1.00 2161 1916
## y1[20] 3.38 3.36 2.18 2.13 -0.26 7.06 1.00 1930 2011
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
ex2-2.stan
fit to poisson dist.
data{
int N;
array[N] int y;
}
parameters{
real<lower=0> l;
}
model{
y~poisson(l);
}
generated quantities{
array[N] int y1;
for(i in 1:N)
y1[i]=poisson_rng(l);
}mdl=cmdstan_model('./ex2-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -4.6355
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 15.2167 8.35723e-05 1.05621e-05 1 1 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 15.22
## l 3.40
## y1[1] 4.00
## y1[2] 6.00
## y1[3] 2.00
## y1[4] 4.00
## y1[5] 3.00
## y1[6] 8.00
## y1[7] 3.00
## y1[8] 4.00
## y1[9] 4.00
## y1[10] 3.00
## y1[11] 3.00
## y1[12] 3.00
## y1[13] 6.00
## y1[14] 0.00
## y1[15] 5.00
## y1[16] 2.00
## y1[17] 4.00
## y1[18] 2.00
## y1[19] 6.00
## y1[20] 2.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 15.96 16.24 0.72 0.29 14.56 16.45 1.00 852 1068
## l 3.43 3.42 0.41 0.39 2.79 4.15 1.00 731 939
## y1[1] 3.52 3.00 1.92 1.48 1.00 7.00 1.00 1689 1711
## y1[2] 3.44 3.00 1.86 1.48 1.00 7.00 1.00 1881 1926
## y1[3] 3.45 3.00 1.86 1.48 1.00 7.00 1.00 1740 1832
## y1[4] 3.40 3.00 1.91 1.48 1.00 7.00 1.00 1919 1988
## y1[5] 3.48 3.00 1.88 1.48 1.00 7.00 1.00 1692 1853
## y1[6] 3.44 3.00 1.91 1.48 1.00 7.00 1.00 2017 1814
## y1[7] 3.50 3.00 1.91 1.48 1.00 7.00 1.00 1716 2033
## y1[8] 3.38 3.00 1.89 1.48 1.00 7.00 1.00 1744 1959
## y1[9] 3.47 3.00 1.93 1.48 1.00 7.00 1.00 1927 2035
## y1[10] 3.38 3.00 1.91 1.48 1.00 7.00 1.00 1710 1892
## y1[11] 3.39 3.00 1.87 1.48 1.00 7.00 1.00 1799 1831
## y1[12] 3.45 3.00 1.91 1.48 1.00 7.00 1.00 1892 1863
## y1[13] 3.42 3.00 1.89 1.48 1.00 7.00 1.00 1838 1913
## y1[14] 3.44 3.00 1.90 1.48 1.00 7.00 1.00 1816 1846
## y1[15] 3.41 3.00 1.92 1.48 1.00 7.00 1.00 1684 1636
## y1[16] 3.33 3.00 1.89 1.48 1.00 7.00 1.00 1805 1880
## y1[17] 3.42 3.00 1.88 1.48 1.00 7.00 1.00 1893 1826
## y1[18] 3.35 3.00 1.84 1.48 1.00 7.00 1.00 1740 1905
## y1[19] 3.46 3.00 1.85 1.48 1.00 7.00 1.00 1868 1917
## y1[20] 3.40 3.00 1.90 1.48 1.00 7.00 1.00 1633 2023
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
difference test
ex3.stan
mean difference of 2 normal distributions
data{
int N;
vector[N] a;
vector[N] b;
}
parameters{
real ma;
real<lower=0> sa;
real mb;
real<lower=0> sb;
}
model{
a~normal(ma,sa);
b~normal(mb,sb);
}
generated quantities{
real d;
d=mb-ma;
}n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)
mdl=cmdstan_model('./ex3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -1244.97
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 26 -24.8565 0.000370989 0.000220677 1 1 43
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -24.86
## ma 9.89
## sa 0.69
## mb 10.57
## sb 1.84
## d 0.68mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -26.71 -26.36 1.51 1.35 -29.77 -24.93 1.01 796 1389
## ma 9.89 9.89 0.18 0.17 9.60 10.19 1.00 2099 1257
## sa 0.76 0.74 0.14 0.13 0.57 1.01 1.00 2129 1503
## mb 10.56 10.56 0.45 0.46 9.81 11.29 1.01 644 596
## sb 2.03 1.99 0.35 0.33 1.54 2.71 1.00 2386 1454
## d 0.66 0.65 0.49 0.49 -0.14 1.50 1.01 727 746
d=mcmc$draws('d')
d## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
##
## chain
## iteration 1 2 3 4
## 1 0.42 -0.15 0.81 0.67
## 2 0.87 1.48 1.15 0.54
## 3 -0.11 1.04 1.30 1.22
## 4 0.89 0.66 1.50 0.58
## 5 0.79 0.98 1.34 0.49
##
## # ... with 495 more iterationsmean(d>0)## [1] 0.911normal regression
ex4-1.stan
single normal regression
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)
data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -1.72607e+06
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 42 -52.917 0.00706613 0.000642451 0.9445 0.09445 117
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -52.92
## b0 14.59
## b1 3.00
## s 8.55mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -52.44 -52.09 1.37 1.16 -55.13 -50.96 1.01 320 508
## b0 14.46 14.62 4.57 3.94 7.08 22.11 1.01 219 167
## b1 3.00 3.00 0.08 0.07 2.87 3.13 1.01 301 398
## s 9.75 9.53 1.74 1.64 7.40 13.02 1.00 901 892
b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)
ex4-2.stan
single normal regression, prediction from new explanatory variable
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N1] m;
vector[N1] y1;
for(i in 1:N1){
m[i]=b0+b1*x1[i];
y1[i]=normal_rng(m[i],s);
}
}# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex4-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -1.88247e+06
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 48 -52.917 0.000245332 0.00357952 0.4582 0.4582 107
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -52.92
## b0 14.59
## b1 3.00
## s 8.55
## m[1] 28.46
## m[2] 58.37
## m[3] 88.28
## m[4] 118.19
## m[5] 148.10
## m[6] 178.01
## m[7] 207.93
## m[8] 237.84
## m[9] 267.75
## m[10] 297.66
## y1[1] 37.39
## y1[2] 61.60
## y1[3] 85.19
## y1[4] 125.25
## y1[5] 150.33
## y1[6] 168.23
## y1[7] 204.06
## y1[8] 236.88
## y1[9] 285.83
## y1[10] 297.70mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.3 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -52.32 -52.00 1.27 0.99 -54.78 -50.94 1.00 580 942
## b0 14.89 14.65 4.15 4.17 8.57 21.85 1.01 226 267
## b1 2.99 3.00 0.07 0.07 2.87 3.11 1.01 303 494
## s 9.71 9.45 1.72 1.58 7.38 12.86 1.00 1038 955
## m[1] 28.75 28.48 3.86 3.86 22.80 35.18 1.01 226 278
## m[2] 58.63 58.43 3.27 3.31 53.45 64.03 1.01 232 303
## m[3] 88.52 88.42 2.75 2.69 84.10 93.07 1.01 253 367
## m[4] 118.40 118.35 2.36 2.30 114.58 122.34 1.01 317 536
## m[5] 148.29 148.29 2.16 2.07 144.70 151.76 1.00 532 879
## m[6] 178.17 178.21 2.21 2.10 174.57 181.63 1.00 1326 1338
## m[7] 208.06 208.04 2.50 2.37 204.00 212.02 1.00 2336 1609
## m[8] 237.94 237.94 2.95 2.75 233.29 242.69 1.00 1851 1428
## m[9] 267.83 267.79 3.50 3.22 262.33 273.36 1.00 1162 1316
## m[10] 297.72 297.70 4.11 3.75 291.17 304.27 1.00 855 1183
## y1[1] 28.45 28.26 10.42 9.93 11.31 46.14 1.00 1016 1737
## y1[2] 58.46 58.79 10.41 10.21 41.27 75.00 1.00 1099 1712
## y1[3] 88.45 88.56 10.27 9.93 71.27 104.58 1.00 1455 1766
## y1[4] 118.69 118.58 9.93 9.56 102.81 135.43 1.00 1561 1924
## y1[5] 148.46 148.51 9.81 9.52 132.63 164.67 1.00 1816 1728
## y1[6] 178.02 178.02 10.07 9.49 161.34 194.70 1.00 1993 2016
## y1[7] 207.81 207.74 10.06 9.85 191.13 224.36 1.00 2058 1725
## y1[8] 237.68 237.82 10.40 10.19 220.94 255.18 1.00 2166 2010
## y1[9] 268.08 268.11 10.63 10.46 250.08 285.64 1.00 1634 2101
## y1[10] 297.88 298.02 10.70 10.19 280.66 315.03 1.00 1519 1784
m=mcmc$draws('m')
summary(m)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 m[1] 28.7 28.5 3.86 3.86 22.8 35.2 1.01 227. 279.
## 2 m[2] 58.6 58.4 3.27 3.31 53.5 64.0 1.01 233. 304.
## 3 m[3] 88.5 88.4 2.75 2.69 84.1 93.1 1.01 253. 367.
## 4 m[4] 118. 118. 2.36 2.30 115. 122. 1.01 318. 536.
## 5 m[5] 148. 148. 2.16 2.07 145. 152. 1.00 533. 879.
## 6 m[6] 178. 178. 2.21 2.10 175. 182. 0.999 1327. 1338.
## 7 m[7] 208. 208. 2.50 2.37 204. 212. 1.00 2336. 1609.
## 8 m[8] 238. 238. 2.95 2.75 233. 243. 1.00 1851. 1429.
## 9 m[9] 268. 268. 3.50 3.22 262. 273. 1.00 1162. 1317.
## 10 m[10] 298. 298. 4.11 3.75 291. 304. 1.00 855. 1184.smm=summary(m)
y1=mcmc$draws('y1')
summary(y1)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 y1[1] 28.5 28.3 10.4 9.93 11.3 46.1 1.00 1017. 1738.
## 2 y1[2] 58.5 58.8 10.4 10.2 41.3 75.0 1.00 1099. 1712.
## 3 y1[3] 88.5 88.6 10.3 9.93 71.3 105. 1.00 1456. 1766.
## 4 y1[4] 119. 119. 9.93 9.56 103. 135. 1.00 1562. 1924.
## 5 y1[5] 148. 149. 9.81 9.52 133. 165. 1.00 1816. 1728.
## 6 y1[6] 178. 178. 10.1 9.49 161. 195. 1.00 1994. 2017.
## 7 y1[7] 208. 208. 10.1 9.85 191. 224. 1.00 2059. 1725.
## 8 y1[8] 238. 238. 10.4 10.2 221. 255. 1.00 2166. 2011.
## 9 y1[9] 268. 268. 10.6 10.5 250. 286. 0.999 1635. 2102.
## 10 y1[10] 298. 298. 10.7 10.2 281. 315. 1.00 1519. 1785.smy=summary(y1)
xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)
par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=m),xy,col='black')+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')
ex4-3.stan
single normal regression, prediction from original explanatory variable
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m;
vector[N] y1;
for(i in 1:N){
m[i]=b0+b1*x[i];
y1[i]=normal_rng(m[i],s);
}
}data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-3.stan')
mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -52.25 -51.94 1.20 0.94 -54.64 -50.95 1.00 509 773
## b0 14.77 14.68 4.03 4.09 8.86 21.52 1.01 177 193
## b1 2.99 2.99 0.07 0.07 2.87 3.10 1.01 244 381
## s 9.61 9.39 1.70 1.60 7.30 12.79 1.00 1039 1025
## m[1] 51.16 51.06 3.35 3.42 46.28 56.72 1.01 182 192
## m[2] 267.85 267.91 3.21 3.12 262.49 273.13 1.00 1377 1216
## m[3] 137.28 137.24 2.21 2.17 133.63 140.89 1.00 337 676
## m[4] 227.00 227.02 2.58 2.52 222.79 231.27 1.00 2040 1404
## m[5] 275.20 275.27 3.34 3.21 269.59 280.67 1.00 1251 1229
## m[6] 74.78 74.65 2.95 3.06 70.39 79.58 1.01 192 205
## m[7] 261.01 261.07 3.10 2.99 255.83 266.09 1.00 1518 1294
## m[8] 28.63 28.51 3.77 3.86 23.19 34.92 1.01 178 194
## m[9] 169.78 169.78 2.14 2.05 166.37 173.32 1.00 773 1077
## m[10] 189.54 189.57 2.21 2.19 185.95 193.19 1.00 1493 1417
## m[11] 48.19 48.08 3.41 3.47 43.24 53.84 1.01 181 199
## m[12] 290.13 290.19 3.61 3.43 284.05 296.04 1.00 1051 1117
## m[13] 181.12 181.13 2.17 2.11 177.66 184.72 1.00 1098 1372
## m[14] 297.48 297.56 3.75 3.57 291.21 303.56 1.00 975 1067
## m[15] 143.22 143.19 2.18 2.13 139.67 146.80 1.00 383 708
## m[16] 160.10 160.09 2.13 2.07 156.62 163.61 1.00 592 1039
## m[17] 65.63 65.51 3.10 3.21 61.05 70.70 1.01 187 190
## m[18] 224.37 224.37 2.55 2.51 220.26 228.61 1.00 2057 1363
## m[19] 77.34 77.19 2.91 3.03 72.98 82.02 1.01 194 220
## m[20] 71.42 71.29 3.01 3.10 66.91 76.28 1.01 190 202
## y1[1] 50.94 50.71 10.43 10.27 34.30 68.01 1.00 1178 1556
## y1[2] 267.63 267.83 10.16 9.59 250.63 284.15 1.00 2146 1856
## y1[3] 137.27 137.50 10.22 9.91 120.33 153.66 1.00 1584 1858
## y1[4] 227.16 227.24 10.37 10.20 209.93 244.14 1.00 2000 1772
## y1[5] 274.90 274.96 10.30 10.17 258.01 291.44 1.00 1715 1649
## y1[6] 75.05 74.89 10.17 10.18 58.26 91.83 1.00 1471 1894
## y1[7] 261.09 261.15 9.94 9.35 244.81 277.62 1.00 1921 1956
## y1[8] 28.65 28.41 10.65 10.02 11.43 46.09 1.00 1028 1740
## y1[9] 169.62 169.62 10.08 9.98 153.09 185.73 1.00 1849 1731
## y1[10] 189.47 189.50 9.96 9.42 172.91 205.61 1.00 1742 1931
## y1[11] 48.17 48.31 10.19 9.95 30.87 64.46 1.00 1214 1608
## y1[12] 289.68 289.67 10.03 9.51 273.53 306.54 1.00 1938 1805
## y1[13] 181.18 181.20 10.08 9.32 164.58 197.39 1.00 2188 2001
## y1[14] 297.13 297.36 10.34 9.90 280.00 314.44 1.00 1688 1711
## y1[15] 143.04 143.08 10.12 9.68 127.03 159.71 1.00 1975 2046
## y1[16] 160.37 160.53 10.12 9.51 143.45 176.60 1.00 1966 1857
## y1[17] 65.59 65.74 10.09 9.82 49.02 82.23 1.00 1011 1431
## y1[18] 224.05 223.88 10.30 9.99 207.18 241.06 1.00 1981 1523
## y1[19] 77.26 77.59 10.33 9.60 60.38 93.23 1.00 1203 1778
## y1[20] 71.77 71.56 10.21 10.16 55.40 88.60 1.00 1503 1567
y1=mcmc$draws('y1')
smy=summary(y1)
par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
estimate correlation
ex4-41.stan
use covariance matrix and multinormal distribution
data{
int N;
array[N] vector[2] y;
}
parameters{
vector[2] m;
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
y~multi_normal(m,cv);
}ex4-42.stan
use covariance matrix and wishart distribution
data{
int N;
matrix[2,2] dp;
}
parameters{
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
dp~wishart(N-1,cv);
}ex4-43.stan
use correlation matrix and wishart distribution
data{
int N;
int M;
matrix[M,M] dp;
}
parameters{
vector<lower=0>[M] s;
corr_matrix[M] r;
}
transformed parameters{
cov_matrix[M] cv;
cv=quad_form_diag(r,s);
}
model{
dp~wishart(N-1,cv);
}ex4-44.stan use covariance matrix and multinormal distribution
data{
int N;
int K;
array[N] vector[K] y;
}
parameters{
vector[K] m;
cov_matrix[K] cov;
}
model{
y~multi_normal(m,cov);
}library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)## [,1] [,2]
## [1,] 0.798 0.626
## [2,] 0.626 0.770cor(y)## [,1] [,2]
## [1,] 1.000 0.798
## [2,] 0.798 1.000plot(y)
#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')
data=list(N=n,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -77.4925
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 26 -3.96424 1.35683e-05 0.000643409 0.4565 0.4565 30
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -3.96
## m[1] 0.08
## m[2] -0.34
## s[1] 0.87
## s[2] 0.86
## r 0.80
## cv[1,1] 0.76
## cv[2,1] 0.59
## cv[1,2] 0.59
## cv[2,2] 0.73mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -8.65 -8.31 1.71 1.52 -11.86 -6.54 1.00 661 864
## m[1] 0.09 0.09 0.23 0.22 -0.29 0.44 1.00 1057 908
## m[2] -0.33 -0.33 0.22 0.21 -0.69 0.02 1.00 985 1080
## s[1] 0.96 0.94 0.16 0.15 0.73 1.25 1.00 1002 926
## s[2] 0.94 0.92 0.17 0.15 0.71 1.25 1.00 957 951
## r 0.76 0.77 0.11 0.09 0.55 0.89 1.01 475 577
## cv[1,1] 0.94 0.88 0.34 0.28 0.53 1.57 1.00 1002 926
## cv[2,1] 0.71 0.65 0.30 0.24 0.34 1.26 1.00 600 647
## cv[1,2] 0.71 0.65 0.30 0.24 0.34 1.26 1.00 600 647
## cv[2,2] 0.91 0.84 0.35 0.27 0.51 1.57 1.00 957 951
#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n
mdl=cmdstan_model('./ex4-42.stan')
data=list(N=n,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -62.2028
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -4.7406 3.36999e-05 4.25106e-05 1 1 21
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -4.74
## s[1] 0.89
## s[2] 0.88
## r 0.80
## cv[1,1] 0.80
## cv[2,1] 0.63
## cv[1,2] 0.63
## cv[2,2] 0.77mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -8.20 -7.85 1.27 1.05 -10.71 -6.84 1.00 708 983
## s[1] 0.95 0.93 0.17 0.15 0.72 1.28 1.00 806 862
## s[2] 0.93 0.91 0.17 0.15 0.71 1.25 1.00 825 1017
## r 0.75 0.77 0.11 0.10 0.54 0.89 1.01 406 475
## cv[1,1] 0.93 0.86 0.34 0.27 0.53 1.64 1.00 806 862
## cv[2,1] 0.69 0.63 0.30 0.24 0.34 1.29 1.00 529 607
## cv[1,2] 0.69 0.63 0.30 0.24 0.34 1.29 1.00 529 607
## cv[2,2] 0.90 0.82 0.34 0.27 0.51 1.56 1.00 825 1017
#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')
m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -36.6064
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -4.7406 0.000412103 1.97173e-05 1 1 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -4.74
## s[1] 0.89
## s[2] 0.88
## r[1,1] 1.00
## r[2,1] 0.80
## r[1,2] 0.80
## r[2,2] 1.00
## cv[1,1] 0.80
## cv[2,1] 0.63
## cv[1,2] 0.63
## cv[2,2] 0.77mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.51 -7.17 1.27 1.05 -10.01 -6.15 1.01 636 858
## s[1] 0.95 0.93 0.17 0.15 0.72 1.26 1.01 1065 934
## s[2] 0.93 0.91 0.17 0.15 0.72 1.24 1.00 1040 715
## r[1,1] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## r[2,1] 0.75 0.77 0.11 0.10 0.54 0.89 1.01 793 530
## r[1,2] 0.75 0.77 0.11 0.10 0.54 0.89 1.01 793 530
## r[2,2] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## cv[1,1] 0.93 0.87 0.34 0.28 0.52 1.59 1.01 1065 934
## cv[2,1] 0.69 0.63 0.31 0.24 0.33 1.26 1.00 805 760
## cv[1,2] 0.69 0.63 0.31 0.24 0.33 1.26 1.00 805 760
## cv[2,2] 0.90 0.83 0.35 0.26 0.52 1.54 1.00 1040 715
mdl=cmdstan_model('./ex4-44.stan')
k=2
data=list(N=n,K=k,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -261.095
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 20 -3.96424 7.74109e-06 0.000432026 0.3788 0.3788 23
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -3.96
## m[1] 0.08
## m[2] -0.34
## cov[1,1] 0.76
## cov[2,1] 0.59
## cov[1,2] 0.59
## cov[2,2] 0.73mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -6.97 -6.60 1.84 1.49 -10.54 -4.79 1.01 706 753
## m[1] 0.07 0.07 0.24 0.23 -0.32 0.46 1.00 926 815
## m[2] -0.34 -0.34 0.24 0.23 -0.71 0.02 1.00 944 779
## cov[1,1] 1.16 1.05 0.51 0.39 0.60 2.04 1.01 1200 952
## cov[2,1] 0.91 0.81 0.44 0.33 0.42 1.70 1.01 1011 717
## cov[1,2] 0.91 0.81 0.44 0.33 0.42 1.70 1.01 1011 717
## cov[2,2] 1.11 1.01 0.48 0.37 0.58 1.94 1.01 1098 711
distributions
normal dist.
distribution of mean, it's from central limit theorem
y~N(m,s), y(-Inf,Inf), E[y]=m, V[y]=s^2ex1.stan
data{
int N;
vector[N] y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~normal(m,s);
}y=rnorm(10,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.91 1.16 1.72 2.04 2.40 4.69par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=length(y),y=y)
mdl=cmdstan_model('./ex1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -10.5958
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 7 -6.44154 4.82669e-05 6.98265e-05 0.8803 0.8803 10
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -6.44
## m 2.04
## s 1.16mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.41 -7.07 1.14 0.81 -9.76 -6.32 1.00 755 813
## m 2.04 2.04 0.48 0.44 1.25 2.83 1.00 752 784
## s 1.43 1.34 0.42 0.34 0.93 2.20 1.00 1514 1324
Bernoulli dist.
The event with probablity p occur y=1, not occur y=0
y~Ber(p), y{0,1}, E[y]=p, V[y]=p(1-p)ex0-0.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
p~beta(1,1);
y~bernoulli(p);
}data=list(N=9,y=sample(0:1,9,replace=T))
mdl=cmdstan_model('./ex0-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -8.44677
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -5.72863 0.000161808 6.74689e-07 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -5.73
## p 0.67mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.70 -7.45 0.68 0.32 -8.96 -7.21 1.00 1035 1383
## p 0.64 0.65 0.13 0.14 0.40 0.84 1.01 880 1045
Poisson dist.
The event occur l times in unit range, the event occur y times.
y~Po(l), y{0,1...Inf}, E[y]=l, V[y]=lex2-2.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0> l;
}
model{
y~poisson(l);
}
generated quantities{
array[N] int y1;
for(i in 1:N)
y1[i]=poisson_rng(l);
}n=10
y=rpois(n,3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 1.25 2.00 2.90 4.50 7.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex2-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -15.7744
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 1.87661 3.45392e-05 2.78157e-06 1 1 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 1.88
## l 2.90
## y1[1] 2.00
## y1[2] 1.00
## y1[3] 2.00
## y1[4] 1.00
## y1[5] 2.00
## y1[6] 2.00
## y1[7] 1.00
## y1[8] 1.00
## y1[9] 4.00
## y1[10] 0.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 2.46 2.72 0.67 0.33 1.20 2.96 1.00 895 1242
## l 2.99 2.96 0.54 0.56 2.15 3.92 1.00 470 830
## y1[1] 2.94 3.00 1.78 1.48 0.00 6.00 1.00 1900 1810
## y1[2] 2.97 3.00 1.76 1.48 0.00 6.00 1.00 1938 1894
## y1[3] 2.97 3.00 1.82 1.48 0.00 6.00 1.00 1964 1794
## y1[4] 3.01 3.00 1.83 1.48 0.00 6.00 1.00 1523 1860
## y1[5] 2.97 3.00 1.81 1.48 0.00 6.00 1.00 1700 1683
## y1[6] 2.94 3.00 1.81 1.48 0.00 6.00 1.00 1538 1798
## y1[7] 3.00 3.00 1.80 1.48 0.00 6.00 1.00 1700 1884
## y1[8] 2.96 3.00 1.80 1.48 0.00 6.00 1.00 1807 1834
## y1[9] 2.99 3.00 1.80 1.48 0.00 6.00 1.00 1732 2021
## y1[10] 2.96 3.00 1.77 1.48 0.00 6.00 1.00 1916 1895
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
binomial dist.
The event with probablity p occur y times on n trials
y~B(n,p), y{0,1...n}, E[y]=np, V[y]=np(1-p)ex2-3.stan
data{
int N;
array[N] int n;
array[N] int y;
}
parameters{
real<lower=0> p;
}
model{
y~binomial(n,p);
}n=10
n0=floor(runif(n,1,10))
y=rbinom(n,n0,0.3)
data=list(N=n,n=n0,y=y)
mdl=cmdstan_model('./ex2-3.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Error evaluating the log probability at the initial value.
## Exception: binomial_lpmf: Probability parameter is 3.32476, but must be in the interval [0, 1] (in '/tmp/RtmpXPrRnk/model-8004714516a.stan', line 12, column 2 to column 18)
## Rejecting initial value:
## Error evaluating the log probability at the initial value.
## Exception: binomial_lpmf: Probability parameter is 1.56971, but must be in the interval [0, 1] (in '/tmp/RtmpXPrRnk/model-8004714516a.stan', line 12, column 2 to column 18)
## Rejecting initial value:
## Error evaluating the log probability at the initial value.
## Exception: binomial_lpmf: Probability parameter is 2.72465, but must be in the interval [0, 1] (in '/tmp/RtmpXPrRnk/model-8004714516a.stan', line 12, column 2 to column 18)
## Initial log joint probability = -27.6751
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -25.4602 0.000303008 0.000379052 0.9528 0.9528 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -25.46
## p 0.36mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -26.97 -26.71 0.70 0.34 -28.36 -26.46 1.00 965 932
## p 0.36 0.36 0.07 0.08 0.25 0.49 1.00 673 767
log normal dist.
logalithm of variable follows normal distribution
log y~N(m,s), y>0ex2-4.stan
data{
int N;
vector[N]<lower=0> y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~lognormal(m,s);
}n=10
y=rlnorm(n,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.9 2.1 4.3 9.8 14.8 34.9par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex2-4.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -37.1816
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -15.534 0.000182535 0.000119359 0.9162 0.9162 20
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -15.53
## m 1.67
## s 1.14mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.50 -16.17 1.10 0.79 -18.67 -15.44 1.00 804 982
## m 1.67 1.65 0.46 0.43 0.96 2.47 1.01 465 625
## s 1.42 1.34 0.41 0.34 0.92 2.17 1.00 1388 1446
exponential dist.
The event occur l times in unit range, the event take time y to occur
y~ex(l), y>0, E[y]=1/l, V[y]=1/l^2ex2-5.stan
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> l;
}
model{
y~exponential(l);
}n=10
y=rexp(n,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.04 0.26 0.92 1.01 1.35 3.19par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex2-5.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -18.0237
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -10.0976 0.00104077 9.54535e-05 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -10.10
## l 0.99mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -10.55 -10.27 0.68 0.29 -12.00 -10.06 1.00 928 1106
## l 1.08 1.05 0.32 0.31 0.61 1.65 1.00 687 873
gamma dist.
The event occur l times in unit range, the event take time y to occur a times
y~Ga(a,l), y>0, E[y]=a/l, V[y]=a/l^2ex2-6.stan
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> a;
real<lower=0> l;
}
model{
y~gamma(a,l);
}n=10
y=rgamma(n,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.11 1.54 2.13 2.28 2.65 5.97par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex2-6.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -25.3383
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -17.9853 0.00100604 1.24513e-06 1 1 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -17.99
## a 1.35
## l 0.59mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -18.98 -18.65 1.06 0.78 -21.17 -17.96 1.01 626 643
## a 1.66 1.57 0.63 0.61 0.80 2.84 1.01 580 785
## l 0.78 0.73 0.33 0.30 0.32 1.40 1.01 549 546
beta dist.
use as prior of binomial dist parameter p
p~Be(a,b), p[0,1], E[p]=a/(a+b), V[p]=ab/(a+b)^2/(a+b+1)ex2-7.stan
data{
int N;
vector[N] p;
}
parameters{
real<lower=0> a;
real<lower=0> b;
}
model{
p~beta(a,b);
}n=20
p=rbeta(n,1,2)
summary(p)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.066 0.163 0.289 0.351 0.441 0.896par(mfrow=c(1,2))
hist(p,main='p')
density(p) |> plot(main='p')
data=list(N=n,p=p)
mdl=cmdstan_model('./ex2-7.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -132.547
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 3.14447 0.000459542 5.50416e-05 1 1 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 3.14
## a 1.24
## b 2.13mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 3.23 3.54 1.01 0.75 1.18 4.21 1.00 1001 1085
## a 1.39 1.35 0.38 0.38 0.83 2.07 1.01 577 997
## b 2.43 2.39 0.69 0.67 1.39 3.64 1.01 604 1025
Dirichlet dist.
use as prior of categorical dist parameter p
p~dir(p0), p[0,1]ex2-8.stan
data {
int N;
int K;
matrix[N, K] p;
}
parameters {
simplex[K] p0;
}
model {
for(i in 1:N){
p[i]~dirichlet(p0);
}
}library(gtools)
n=20
p=rdirichlet(n,c(0.5,0.3,0.2))
summary(p)## V1 V2 V3
## Min. :0.002 Min. :0.000 Min. :0.000
## 1st Qu.:0.120 1st Qu.:0.007 1st Qu.:0.006
## Median :0.547 Median :0.117 Median :0.033
## Mean :0.542 Mean :0.328 Mean :0.130
## 3rd Qu.:0.930 3rd Qu.:0.673 3rd Qu.:0.166
## Max. :0.998 Max. :0.998 Max. :0.896boxplot(p)
data=list(N=n,K=ncol(p),p=p)
mdl=cmdstan_model('./ex2-8.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = 51.2221
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 53.1332 0.00825801 0.00122163 1 1 6
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 53.13
## p0[1] 0.49
## p0[2] 0.27
## p0[3] 0.23mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 48.70 49.01 0.98 0.71 46.73 49.63 1.00 1054 1197
## p0[1] 0.48 0.48 0.06 0.06 0.39 0.58 1.00 1965 1405
## p0[2] 0.28 0.28 0.05 0.05 0.20 0.37 1.00 1412 1167
## p0[3] 0.24 0.24 0.04 0.04 0.17 0.31 1.00 1726 1408
categorical dist. with parameters following Dirichlet dist.
h~Dir(a), h[0,1], sum(h)=1, a[0,1], sum(a)=1
y~Cat(h), y[0,1], sum(y)=1ex2-9.stan
data {
int N;
int K;
array[N] int<lower=1,upper=K> y;
}
parameters {
simplex[K] a;
simplex[K] h;
}
model {
h~dirichlet(a);
y~categorical(h);
}library(gtools)
n=20
h=rdirichlet(n,c(0.5,0.3,0.2))
summary(h)## V1 V2 V3
## Min. :0.015 Min. :0.000 Min. :0.000
## 1st Qu.:0.126 1st Qu.:0.025 1st Qu.:0.004
## Median :0.253 Median :0.214 Median :0.087
## Mean :0.422 Mean :0.297 Mean :0.281
## 3rd Qu.:0.758 3rd Qu.:0.554 3rd Qu.:0.593
## Max. :1.000 Max. :0.974 Max. :0.985c0=c(1,2,3)
for(i in 1:n) y[i]=sample(c0,1,prob=h[i,])
table(y) |> prop.table()## y
## 1 2 3
## 0.45 0.35 0.20data=list(N=n,K=ncol(h),y=y)
mdl=cmdstan_model('./ex2-9.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -28.4628
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 13 -21.589 0.000450354 2.27097e-05 1 1 17
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -21.59
## a[1] 0.37
## a[2] 0.34
## a[3] 0.28
## h[1] 0.47
## h[2] 0.35
## h[3] 0.18mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -30.60 -30.24 1.54 1.38 -33.46 -28.79 1.01 844 1062
## a[1] 0.34 0.33 0.18 0.19 0.08 0.66 1.00 1379 1261
## a[2] 0.34 0.33 0.18 0.19 0.08 0.68 1.00 1213 983
## a[3] 0.31 0.29 0.17 0.18 0.06 0.61 1.00 1061 753
## h[1] 0.45 0.45 0.11 0.11 0.27 0.63 1.00 3059 1629
## h[2] 0.35 0.34 0.10 0.11 0.19 0.53 1.00 2371 1592
## h[3] 0.20 0.19 0.09 0.09 0.08 0.37 1.00 2393 1689
